Towards a theoretical description of molecular junctions in the Coulomb blockade 

regime based on density functional theory 
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Non-equilibrium Greens function techniques (NEGF) combined with Density Functional Theory 
(DFT) calculations have become a standard tool for the description of electron transport through 
single molecule nano-j unctions in the coherent tunneling regime. However, the applicability of these 
methods for transport in the Coulomb blockade (CB) regime is still under debate. We present 
here NEGF-DFT calculations performed on simple model systems in the presence of an effective 
gate potential. The results show that: i) the CB addition energies can be predicted with such 
an approach with reasonable accuracy; ii) neither the magnitude of the Kohn-Sham gap nor the 
lack of a derivative discontinuity in the exchange-correlation functional represent a problem for this 
purpose. 

PACS numbers: 71.10.-w, 71.15.Fv, 71. 15. Mb, 73.23.Hk 
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Interest in electron transport between nanoscale con- 
tacts has recently intensified, due to the advent of the 
technologically motivated field of molecular electronics 
and recent progress in experimental techniques for ma- 
nipulating and contacting individual molecules^. Elec- 
tron transport is generally described within two limit- 
ing regimes, namely, coherent transport (CT) for strong 
coupling between the molecule and the electrodes or 
Coulomb blockade (CB) for weak coupling. These 
two regimes can be considered as complementary in a 
quantum-mechanical sensed: CT deals with evanescent 
waves passing through the junction and does not lead to 
a stepwise charging of the central molecule; in the CB 
(incoherent) regime, the electrons cross the junction one 
by one and get localized during a finite time over the 
molecule. Transport in the CB regime is characterized 
by so-called CB diamonds depicting frontiers between 
low- and high-conductivity domains in the source/drain 
bias gate voltage coordinates^^. Key quantities in this 
stability diagram are the energy differences between the 
ionization and affinity levels of a quantum dot or single 
molecule, which are referred to as addition energies E^. 

First-principle non-equilibrium Green's function 
(NEGF) methods are typically implemented^^ in 
combination with density functional theory (DFT) to 
theoretically describe electron transport through single 
molecule junctions. This approach has been used in 
numerous works and proves to be very useful for the 
description of the CT regime. On the other hand, for 
electron transfer in the CB regime, an integer charge 
is transferred between the molecule and the leads and 
results in a relaxation of the electronic structure of 
the central molecule. It has been recognized that this 
process cannot be suitably described with DFT or 
spin-restricted Hartree-Fock (RHF) based ground-state 
calculations— and that a multi-determinant configu- 
ration interaction (CI) scheme, or Fock space, should 
provide a general solution to this problem^. However, 
it is difficult to apply such a CI treatment to the open 



system of the whole molecular junction, especially 
because the leads are usually described by mean-field 
band structure calculations. Recently, the CB regime 
was described within a NEGF approach in conjunction 
with spin-polarized hybrid DFT 11 ^ and Hartree-Fock^ 
methods for junctions based on finite clusters containing 
partially filled degenerate orbitals. Recent quasiparticle 
calculations based on the GW approximation were found 
to not accurately represent the impact of local spin and 
charge fluctuations on CB^i. 

It has also been suggested that a standard DFT frame- 
work has inherent problems for describing electron trans- 
port in both the CB and CT regime^, due to the self- 
interaction (SI) of electrons^ in a Kohn-Sham frame- 
work and the lack of a derivative discontinuity (DD)i£ in 
the evolution of the eigenenergy of the highest occupied 
molecular orbital (HOMO) as a function of its occupancy, 
which can be fractional. In Ref. 1 -, these two issues have 
been portrayed as intimately linked and a SI correction 
scheme has been devised as a remedy. Independently, 
optimized effective potentials and hybrid methods have 
been employed for reducing or removing SI from DFT- 
based electron transport calculations^. 

In this work, we present NEGF-DFT calculations per- 
formed in the presence of an effective gate potential 
V gate, which allows for an explicit charge transfer be- 
tween one-dimensional lithium chains as electrodes, and 
H2 and benzene molecules as the central unit (see Fig. [I}. 
The results obtained for these model systems show that: 

i) the size of the energy gap estimated from the elec- 
tronic eigenenergies of the frontier Kohn Sham (KS) or- 
bitals is not directly related to the threshold voltages or 
addition energies E af jd associated to the CB diamonds; 

ii) quantitatively realistic values for Fi a dd can be derived 
despite the fact that our calculations are carried out in 
the local density approximation (LDA) for the exchange- 
correlation functional (XC) and introduce a distortion 
in the charging process characteristic for DFT (in the 
sense that distinct one-electron transfer steps are re- 
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FIG. 1: (Color online) Geometry and pattern of the applied 
gate potential in the two model molecular junctions of our 
study: a H2 (top) and a benzene molecule (bottom) are at- 
tached to one-dimensional lithium electrodes, with a fixed 
separation of 4.37 and 4.00 A, respectively. The profile of 
the applied gate potential is shown with grey shading in both 
junctions, as obtained from the difference in the electrostatic 
potential between calculations performed with a V gate of IV 
and 0V; the maximum of the potential is located in the black 
regions. 



placed by continuous curves due the introduction of frac- 
tional charges N add ). 

Fig. Q] illustrates the two systems of our study, for 
which we performed NEGF-DFT calculations with the 
commercially available ATK software^. A single Li atom 
is included in the scattering region on each side of the 
central unit and each electrode part contains six atoms. 
A self-consistent solution for the electron density of the 
open system as a whole is carried out within the Keldysh 
formalism^iiiSii^ for every value of V gate on a grid of 0.5 
V; a small source-drain bias of 20 mV has been applied 
in order to generate a current through the junction. The 
voltage V 'gate is introduced in the Keldysh Hamiltonian 
as H^ v = H^ v + VgateSfw, where is the overlap ma- 
trix and the indices fi and v run only over the atomic 
basis functions of the molecules (i.e. N gate is not applied 
to the electrodes); the term containing V gate is added 
to the Hamiltonian at every step in the self-consistent 
cycle. The shape of the resulting potential is shown as 
gray shades in Fig. [TJ In all calculations presented in 
this article LDA has been used for the XC functional 
and a double-zeta polarized (DZP) basis set for both the 
molecules and the electrodes. 

Fig. [5] shows the evolution of the charge on the 
molecules (N add ) and the KS-MO eigenenergies as a func- 
tion of V gate for both systems. The range of gate voltages 
is larger in the positive region than in the negative one 
due to the well-known convergence problems for anions 
when using localized basis sets^i. In contrast to the cor- 
rect physical behaviour of electrons moving between the 
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FIG. 2: (Color online) (Top) Evolution of the number of elec- 
trons on the molecules (N a dd) as a function of V 'gate together 
with a schematic description of the CB diamonds at the pre- 
dicted threshold voltages for the junction with H2 (left) and 
benzene (right), respectively; (Bottom) The corresponding 
evolution of the MO eigenenergies as a function of Vgate. 
N a dd has been determined by integrating charge density dif- 
ferences, with V gate = 0V taken as a reference and the spatial 
border between the molecule and the electrode defined as the 
region invariant in charge to V gate. The MO eigenenergies 
have been extracted from NEGF-DFT calculations by using 
a projection scheme— and by defining the zero energy at the 
Fermi level of the Li electrodes. 



electrodes and the molecule one by one at given threshold 
voltages, N add evolves linearly from to almost 2 elec- 
trons with Vgate due to the use of a closed-shell ansatz 
and the lack of a derivative discontinuity in DFT— 
This general deficiency of DFT allows for any fractional 
value between and 2 for the occupation of frontier or- 
bitals of the molecules, even if they are only weakly cou- 
pled to metal electrodes (as it is the case in the present 
calculations). 

Nevertheless, the addition energies Eadd associated to 
the CB diamonds can be calculated as: 

E add = (E(N + 1) - E(N)) - (E(N) - E(N - 1)) = 

/ N add Vgate(N add ) - / N add V gate(N add ) = 
Jo J-l 

Vgate(N add = +0.5e) - Vgate(N add 



0.5e) (1) 



In Eq. [TJ the first line is a general expression to calcu- 
late the difference between the redox levels from the total 
energies of the neutral and charged isolated molecules (as 
we did hereafter using Gaussian™). In a second step 
threshold voltages obtained from NEGF-DFT calcula- 
tions are introduced, where a midpoint rule for the in- 
tegration over the charge has been used^; Fj add is thus 
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TABLE I: Addition energies E ac jd calculated as V gate 
(N add =+0.5e)-Vgate (N add =-0.5e) from NEGF-DFT for the 
junctions shown in Fig. 1, and as E(N+1)+E(N-1)2 E(N) 
from the total energies of the neutral and charged isolated 
molecules. The KS HOMO-LUMO gap extracted from the 
NEGF-DFT calculations and reliable values for Eadd of the 
isolated molecules provided by experiment for flz^S and GW 
calculations for benzene^ are also given for comparison. All 
values are given in eV. 



Molecule 


NEGF-DFT 


BFT-Etotal 


KS-MOs 


Lit.(IP-EA) 


H 2 


21.30 


21.35 


11.95 


18.56 


benzene 


9.80 


11.54 


5.23 


10.51 



evaluated from the gate voltages required to add or sub- 
tract half an electron to the molecule. The second part 
of Eq. [1] is general in electrostatics and makes use of the 
relationship between a total energy E, a voltage V and a 
charge N2£. 

Table Q] collects the values of E ae &j calculated for both 
junctions using the two expressions in Eq. [TJ with the 
same basis set and XC functional. For comparison, we 
also provide the Kohn-Sham HOMO-LUMO gap which 
has been derived from NEGF-DFT calculations at Vgate 
= V by projecting the eigenstates of the semi-infinite 
junction on the sub-Hamiltonian defined by the part 
of the basis set that is localized exclusively on the 
molecule^; reliable values for E a dd are also listed for the 
isolated H2 26 and benzene 2 ^ molecules. N a dd has been es- 
timated by integrating the voltage-dependent differences 
in the self-consistent charge density over the molecular 
part (see the caption of Fig. [2] for more details). The 
counter charges are located on the Li chains, with the 
main part on the atoms directly in contact with the 
molecule and the remainder decaying into the wires in 
Friedel oscillations (not shown here). We stress that the 
scaling factor a used in experiments^ to distinguish be- 
tween the applied gate voltage and the effective potential 
felt by the molecule is equal to one in our calculations. 
Although the KS-MO gap is as expected about a factor 
of two smaller than all the other values, Table [T] shows 
that the results obtained from the NEGF-DFT approach 
match the two reliable sets of "E a dd values obtained for 
the isolated molecules rather well. 

The evolution of the MO eigenenergies as a function of 
the gate voltage in Fig. [2] reflect the corresponding curves 
obtained for N a dd since the occupation of a given MO 
can only vary between and 2 while it stays pinned to 
the Fermi level. ~N a dd thus changes with Vgate while the 
eigenenergies stay constant and vice versa. The pinning 
of the MOs to Fermi level EF is also reflected in the trans- 
mission functions shown in Fig. [3] for the LUMO peak of 
the Li-H2-Li junction. As long as the energy of the trans- 
mission peak does not reach Ep, it shifts as a function 
of Vgate with a slope of one; this evolution is stopped 
when progressive charging is initiated in our mean-field 
calculations, instead of the correct stepwise charging, as 
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FIG. 3: (Color online) Transmission functions for the Li-tb-Li 
junction for different values of Vgate, where the Fermi level 
Ef of the Li chain defines zero on the energy axis. The peak 
corresponding to the LUMO moves from ~+5 eV at Vgate = 
0V to ~+0.2 eV at Vgate = -5V; for Vgate = -10V, it becomes 
pinned to Ef- In the transmission functions obtained for the 
two finite values of Vgate, an additional peak associated to a 
higher-lying unoccupied MO can also be seen above 6 eV. The 
inset shows the current in the same junction with an applied 
source-drain bias of 20 mV as a function of Vgate. 

already discussed above. As a consequence, the peak in 
the current / Vgate curve shown in the inset of Fig. [3J 
is unphysically broadened by up to 10 V; in reality, this 
peak is expected to be much narrower due to the weak 
coupling between the molecule and the electrodes. The 
inset does show the high- and low-conductivity domains 
characteristic of the CB behavior when the gate voltage 
is modulated; the point is, however, that the position of 
the boundaries along the Vgate axis is distorted by the 
mean-field character of the calculations and should be 
extracted properly in the way illustrated in Fig. [2] above. 

It would be highly desirable to conduct similar calcula- 
tions for more realistic electrodes than Li chains and for 
junctions where the molecule is chemisorbed on the leads 
via anchoring groups and the weakness of the coupling 
required for the CB regime promoted by alkyl spacers^. 
There are two main reasons preventing us to perform 
such calculations: i) the rather crude way used to intro- 
duce Vgate by shifting the Keldysh Hamiltonian matrix 
elements in a defined region does not allow for the in- 
troduction of chemical bridges without ambiguity; and 
ii) due to the rather high voltages required for molecules 
with a moderate gap size and in view of the resulting 
significant perturbation of the electronic structure of the 
systems, the convergence of our calculations can only be 
achieved in a stepwise fashion, i.e. by using the con- 
verged density of a given voltage as a starting point for a 
next iteration, with steps of 0.5 V. Since this ultimately 
implies a large number of single self-consistent field cal- 
culations, replacing one-dimensional metallic chains by 
realistic surface electrodes is beyond our current com- 
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putational resources. In this context, we cannot at this 
stage address the issue of screening or image charge ef- 
fects which are highly relevant for the interpretation of 
recent experimental data. This is because such effects 
are expected to be only significant when considering large 
surface areas instead of surfaces consisting of just a sin- 
gle atom on each side. A recent comparison between 
DFT-LDA versus GW results for the adsorption of ben- 
zene on graphite^! suggests that DFT is likely to fail in 
describing such polarization effects. We note, however, 
that the adsorbed molecule was not explicitly charged in 
this study and the focus was on the KS-MO gap, which 
is fundamentally different from our approach. 

In conclusion, we have addressed the issue whether 
NEGF-DFT techniques, commonly used for the theoreti- 
cal description of electron transport in the coherent tun- 
neling regime, can also be employed for deriving relevant 
parameters in the Coulomb blockade regime. In our ap- 
proach, we introduced explicitly a gate voltage in order to 



induce a charge transfer from the electrodes (Li chains) 
to the inserted molecules (H 2 and benzene). We obtained 
realistic results for addition energies (a key quantity in 
CB experiments) by using NEGF-DFT for the calcula- 
tions on these simple model systems; we also argued that 
our results are not affected by the widely known short- 
comings of DFT such as the lack of a derivative discon- 
tinuity in the exchange-correlation functional. 
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